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In recent years many procedures have been proposed to check the anisotropy 
of a dataset. We present a new simple procedure, based on a scale dependent 
approach, to detect anisotropy signatures in a given distribution with particular 
attention to small dataset. The method provides a good discrimination power 
for both large and small dataaets, even in presence of strong contaminating 
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isotropic background. We present some applications to simulated datasets of 
events to investigate statistical features of the method and present and inspect 
its behavior under both the null or the alternative hypothesis. 

Keywords: Anisotropy, Multiscale analysis, Cosmic Rays 

1. Introduction 

In many field involving data analysis, the search for anisotropy has played 
a crucial role. Many experimental data such as cosmological, astrophysical, 
atmospherical or geophysical ones require the search for clustering of objects 
and the measure of deviation from isotropy of a given angular distribution. 
During the last decades, many estimators, namely correlation function [1-3], 
have been proposed. These methods apply to angular coordinates of objects 
as well to distributions of arrival directions of events. 

Different estimators have been defined in the common anisotropy analy- 
sis, by Peebles, Davis-Peebles, Landy-Szalay and Hamilton [1-3]. Recently, 
new estimators have been introduced to study the anisotropy signature of 
a given direction distributions: the modified two-point Rayleigh [4], and 
shape-strength method derived from a principal component analysis of 
triplets of events [5] . 

Within the present work, we present a new fast and simple method for 
anisotropy analysis, which makes use of a multiscale approach and depends 
on one parameter only. The main advantage of our estimator is the possibil- 
ity to analytically treat the results: the analytical approach drastically re- 
duces computation time and makes available the possibility of applications 
to very large datasets. We test the method on several simulated isotropic 
and anisotropic arrival direction distributions and perform an extensive 
analysis of its statistical features under both the null and the alternative 
hypotheses. 

2. MAF: Multiscale Autocorrelation Function 

Let »S be a region of a spherical surface defining a sky and let Pi (</>, 6) {i = 
1, 2, n) be a set of n points on S. The sky S is partitioned within a grid 
of N equal-area (and almost-equal shape) disjoint boxes Bk {k = 1, 2, A^) 
as described in Ref . [6] . Each box Bk covers the solid angle 



where 2Q is the apex angle of a cone covering the same solid angle. Let J7 
be the solid angle covered by 5: N,<d and are deeply related quantities 
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that define a scale. 

The function A{@) that quantifies the deviation of data from an 
isotropic distribution at the scale 9, is chosen to be the Kullback-Leibler 
divergence [7,8] 



where ■i/'/c(0) is the fraction of points in the dataset that fall into the box 
Bk and V'fcC©), generally a function of the domain meshing, is the expected 
number of points on Bk from an isotropic distribution of n points on S. 
The Kullback-Leibler divergence is an information theoretic measure that 
quantifies the error in selecting the density tp{&) to approximate the density 
ip{&) and it is strictly connected to maximum likelihood estimation (see 
Appendix A in Rcf. [13]). 

It is straightforward to show that A{Q) is minimum for an isotropic 
distribution of points or, in general, when ip{Q) ~ 4>{0). If ^data(9) and 
Ai^o(0) refer to the data and to an isotropic realization with the same 
number of points respectively, we define multiscale autocorrelation function 
(MAF) the estimator 



where (Aiso(6)) and (7Ai^^{&) are the sample mean and the sample standard 
deviation, estimated from several isotropic realizations of the data. If 7io 

denotes the null hypothesis of an underlying isotropic distribution for the 
data, the chance probability at the angular scale &, properly penalized 
because of the scan on 6, is the probability 



obtained from the fraction of null models giving a multiscale autocorrela- 
tion, at any angular scale Q' in the parameter space V, greater or equal 
than that of data at the scale 6*. Under the null hypothesis Ho, the es- 
timator s{£)) follows an half-gaussian distribution, independently on the 
value of the angular scale and on the number of events on S, as it will 
be successively shown in the text. 




(1) 



|^data(9)-(Ai3o(9))| 



(2) 



p(e) = Pr(si,o(e') > sdata(e)|Ho,ve' e v) 



(3) 



''The null hypothesis is eventually rejected in favor of the alternative Hi = —•Ho — being 
the negation operator — at the angular scale ©, with probability 1 — p(6). 
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(a) Static counting (b) Dynamic counting 



Fig. 1. A cluster of three points falling into different boxes. 

3. Dynamical Counting 

The simplest definition of the counting algorithm involves the fixed grid 
introduced in Ref. [6]. In some cases, such a static counting approach could 
not reveal an existing cluster. For instance, Figure 1(a) shows a typical 
scenario where some points of a given triplet fall into different cells. Indeed, 
the fixed grid may cut a cluster of points within one or more edges. 

The smoothing, adopted in our study, deals with a new counting proce- 
dure for the estimation of the density ■0^(0). 

From the starting distribution of points, we define a new distribution 
defined as follows: for each point Pi, we consider a set of 8 points distributed 
around Pi as described in Figure 1(b), being the distance between each of 
this points and Pi equal to ^- Within this framework, from each point Pi 
we define a set of 9 points {Pi and the new 8 points), and assign to each 
point the weight i. Finally, we follow the standard procedure described in 
Section 2 by using the weighted distribution of points'^. 

Our numerical studies show that such a dynamical counting approach 
gives more information on the amount of clustering in the data. 

The main difference between the static and the dynamical counting 
lies in the value of the estimator when the procedure is applied to Monte 
Carlo realizations of the sky. For instance, let us consider the Figure 2, 
where we show a clustered (Figure 2(a)) and an unclustered (Figure 2(b)) 
set of points. The static counting is not able to recover the differences 
between the two configurations. Conversely, if the dynamical counting is 
applied, the extended points in Figure 2(a) are concentrated mainly in two 



''In general, in the weighting procedure, we have to take into account many factors, as 
for example the exposure of the experiment. See Ref. [13] for more details. 
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(a) Clustered events (b) Unclustered events 

Fig. 2. a) Three clustered points: the weighted points are concentrated in two adjacent 
boxes, b) Three unclustered points: the weighted points are distributed on the neighbor 
cells. 

adjacent boxes while in Figure 2(b) they are distributed on the neighbor 
cells, this lead to two different Monte Carlo skies producing the same 

clustered configuration shown in Figure 2(a), and of consequence the same 
weight distribution, are not frequently expected: in this case, the value of 
s{Q) should be greater than that one estimated from the static method. 
The direct consequence of a greater value of the estimator s[Q) is a lower 
chance probability and the main advantage of using the dynamical counting, 
instead of the static one, should be the lowest penalization of s(8) only if 
an anisotropy signal is really present. 

In order to illustrate the importance of dynamical counting in 
the anisotropy signal detection, we have generated 5000 isotropic and 
anisotropic skies of 100 events each. In each anisotropic sky, 60% of events 
are normally distributed, with dispersion p, around 10 random sources and 
40% of events are isotropically distributed. For each angular scale Q, we 
have estimated the average value of s(0) with the static and the dynamical 
counting, separately. Results are shown in Figure 3 for p 5° (a), p = 10° 
(b), p = 20° (c) and for the isotropic map (d). In the case of the anisotropic 
skies, the dynamical counting provides a greater estimation of s{Q) than 
the static counting, leading to a smaller estimation of the corresponding 
chance probability and improving the signal-to-noise ratio. 

4. Interpretation of MAF 

In the specific case of MAF, the angular scale Q* , where the significance 
is minimum, turns to be the significative clustering scale: it is the scale 
at which occurs a greater accumulation of points respect to that one oc- 
curing by chance, with no regard for a particular configuration of points, 
e.g. doublets or triplets. To illustrate better the clustering scale detection 
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Fig. 3. MAF: average s{@) (solid line) estimated from 5000 isotropic and anisotropic 
skies of 100 events each. In each anisotropic sky, 60% of events are normally distributed, 
with dispersion p, around 10 random sources and 40% of events are isotropically dis- 
tributed. The dashed line indicates the value of the dispersion adopted to generate the 
corresponding mock map: a) 5, b) 10 and c) 20 degrees; d) isotropic map. 
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Fig, 4, MAF: average chance probability (solid line), with 68% region around the mean 
value, estimated from isotropic and anisotropic skies generated as explained in Figure 
3, Dynamical counting is used. The dashed line indicates the value of the dispersion 
adopted to generate the corresponding mock map: a) 5, b) 10 and c) 20 degrees; d) 
isotropic map. 



feature of MAF, we have generated 5000 isotropic and anisotropic skies of 
100 events each, as previously described at the end of Section 3. Figure 4 
shows the average chance probabihty, with 68% region around the mean 
value, versus the angular scale for three values of the dispersion, namely 
p = 5° {a.), p — 10° (b), p = 20° (c), and for the isotropic map (d). As ex- 
pected, chance probability is close to one and nearly flat in the case of the 
isotropic map, because all clustering scales are equally likely. Conversely, 
for all anisotropic maps, the average chance probability gets a minimum 
around the corresponding value of p. Thus, our estimator is able to recover 
the most significant clustering scale. It should be remarked that when the 
20° dispersion is used, the angular scale of the minimum is less obvious be- 
cause of the large fluctuations due to the isotropic contamination. Finally, 
it is worth noticing that we have observed that the curve around the value 
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of p gets narrower by increasing the number of events. 
5. Statistical analysis of MAF 

We have investigated the statistical features of MAF by inspecting its be- 
havior under both the null or the alternative hypothesis (see Ref. [13] for 
further information). In particular, we have estimated the significance a 
(or Type I error), i.e. the probability to wrongly reject the null hypothesis 
when it is actually true, and the power 1 — /3 (where /? is known as Type II 
error), i.e. the probability to accept the alternative hypothesis when it is in 
fact true. In the following we will adopt the dynamical counting previously 
discussed. 

Null hypothesis. Wc generate isotropic maps of 10^ skies, by varying 
the number of events from 20 to 500 and estimate the MAF for several 
values of the angular scale 6. By choosing: 

p(6*) = argniinp(6) 

We find an excellent flat distribution of probabilities p(0*), as expected 
for analyses under the null hypothesis 'Hq. In other words, MAF is not 
biased against Hq, as required for good statistical estimators. 

Because of the definition in Eq. (1) and of the central limit theorem, 
a Gaussian distribution is expected for the function ^4(9), and of conse- 
quence, the half-normal distribution 

gi/2[5(e)]= e ^ (4) 

for <j{Q) = 1, is expected for the estimator s(0) defined as in Eq. (2), 
being normalized to zero mean and unitary variance. It follows that the 
(unpenalized) probability to obtain by chance a value of the MAF, greater 
or equal than a given value sq, is just 1 — erf(-^), being erf the standard 
error function, independently of the angular scale 0. 

Although this nice feature of the MAF estimator it is more important 
to identify the distribution of the penalized probability p(0), because of 
the required penalization due to the scan over the parameter Q. For this 
reason, the distribution of Smax = max{s(0)} is of interest for applications. 

Our numerical studies show that such a distribution exists and it cor- 
responds to one of the limiting densities in the extreme value theory (see 
Appendix B in Ref. [13]). 
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Fig. 5. MAF: distributions of max{s(0)} for n = 40,60,80, 100 and 500 events. Solid 
line correspond to the least-square fit of the Gumbel density with parameters /i = 1.743± 
0.002 and a = 0.470 ± 0.002 (x^/ndf = 1.1 X 10"^). 

In Figure 5 are shown the probabihty densities of Smax for n = 
40, 60, 80, 100 and 500 events: independently on n, each density is in ex- 
cellent agreement with the Gumbel distribution of extreme values, for the 
parameters /i = 1.743±0.002 and a — 0.470±0.002. Such values correspond 
to the mean and to the standard deviation of the distribution, fi 2.00 and 
cr « 0.59, respectively. It follow that the probability to obtain a maximum 
value of s(0), at any angular scale 6, greater or equal than a given value 
max{s(9)} is 

p(max{s(e)}) = l-exp[-exp( L^-^ — ^1)], 

a 

Alternative hypothesis. In order to investigate the behavior of MAF 
under the alternative hypothesis, we have generated anisotropic maps of 10^ 
skies, by varying the number of events from 20 to 100. 

In order to estimate the power of MAF, we build reasonable anisotropic 
maps of CR events reflecting in part the real-world scenario, keeping in 
mind that our purpose is to build an anisotropic set of events for statistical 
analysis and not to generate events mimicking real datasets with the best 
available approximation (see Ref. [13] for further information). 

A sky is labelled as anisotropic if, for a fixed value of the significance 
a, the penalized chance probability as defined in Eq. (3) is lesser or equal 
than a, i.e. if the condition 

p(8*) — argminp(0) < a 

holds for some angular scale Q*. In Figure 6 is shown the power for two 
values of the significance threshold, namely a — 0.1% and a — 1%, esti- 
mated through the analytical approach. For applications, a power of 90% is 
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Fig. 6. MAF power vs. the number of events sampled from anisotropic mock maps 
generated as described in the text, for values of the significance corresponding to a = 
0.1% and a = 1%. 

generally required: under this threshold the method could miss to detect an 
existing anisotropy signal. In the case of the MAF, and for the considered 
anisotropic mock map, the power increases with the number of events n 
and it is able to detect the anisotropic signal for n > 60, with significance 
a — 1%. However, by decreasing the significance for the statistical test, the 
power requires a greater number of events to reach the 90% threshold, as 
expected. 

6. Conclusion 

We introduced a new statistical test, based on a multiscale approach, for de- 
tecting an anisotropy signal in the arrival direction distribution of UHECR. 
We showed that our procedure is suitable for the analysis of both small and 
large datasets of events. 

The advantages of this new approach over other methods are multiples. 
First, the method allows an analytical description of quantities involved 
in the estimation of the amount of anisotropy signal in the data, avoiding 
thousands of Monte Carlo realizations needed for the penalizing procedure 
of results. Second, the method allows the detection of a physical observable, 
namely the clustering scale. Third, the method is unbiased against the null 
hypothesis and it provides a high discrimination power even in presence 
of strong contaminating isotropic background, for both small and large 
datasets. 

This method is suitable for the detection of the anisotropy signal in 
each data set involving a distribution of angular coordinates on the sphere, 
and it can be adapted to non-spherical spaces by properly redefining the 
dynamical counting algorithm. 
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